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TECHNICAL PAPER 


THE VARIATION OF CORROSION POTENTIAL WITH TIME FOR 
COATED METAL SURFACES 

INTRODUCTION 


The measurement of corrosion potentials as determined using the EG&G-PARC model 

350A Corrosion Measurement Console are being explored and developed. The purpose is to rank and 
select materials and coatings for use in seawater/salt air environments such as those experienced by the 
solid rocket booster cases of the Space Shuttle. 

The variation of corrosion rates, electrical resistances, and polarization resistances, as determined 
by electrochemical methods, with time have been reported earlier [ 1 ] . These results were obtained in a 
study of the behavior of primer-coated 2219-T87 aluminum. Although brief mention of the variation of 
corrosion potential with time was made, no detailed study of this problem was undertaken. In the 
present work, a study of the corrosion potentials for coated 4130 steel and primer coated 2219-T87 
aluminum has been carried out. The coated steel samples were included to provide a comparison with 
the work of Wormwell and Brasher [2], who investigated the effects of surface preparation and the type 
and number of coats of paint on the potential of coated steel specimens immersed in synthetic seawater. 
Their work extended over a much longer time period (about 5 months), and the measured potentials 
represented the average for a large number of measurements. The objective of the present work was to 
make the measurements over a much shorter period of time (about 1 month), and to use only a single 
sample rather than the average for a large number of samples. However, considerable scatter of the 
observed data, particularly where the corrosion potentials were recorded manually, make it necessary to 
develop a computer program to smooth the data before proper interpretation could be made. Details 
of this program are described in the next section, and a complete listing of the FORTRAN-77 Program 
is included in the appendix. 


THE SMOOTHING PROGRAM 


The basis for the data smoothing procedure has been described by Lanczos [3], and involves 
smoothing in the large by Fourier analysis. The entire set of data is treated as one unified whole. The 
routine serves as a low pass filter, in which the true course of the function and superimposed noise are 
separated. This method of smoothing has the advantage that it is more independent of any special 
assumptions concerning the nature of the unknown /(x). Only essential details of the method will be 
given here, and the reader is referred to the original report for further information. 

A large number of observations at the points 


x = 0, h, 2h, ..., nh = £ 


( 1 ) 


is to be considered. If a properly chosen a+|3x is subtracted from /(x), 



g(x) = /(x) - (o+/3x) 


( 2 ) 


where a and j3 are determined by the boundary conditions, 


g(0) = 0, g (2) = 0 


(3) 


and 


g(-x) = -g(x) 


(4) 


the result is a function which, if made periodic with the period 22, has no discontinuity in either func- 
tion or derivative. The function g(x) is developed into a pure sine series of the form: 

g(x) = bj sin j x + b 2 sin-^p- x + ... . ( 5 ) 


Since 


y k = /(kh), (k = 0, 1, 2, ..., n) 


( 6 ) 


/(x) is modified to 


g(x) = /(x) - f(0) - 


m ~ /(0)x 

2 


(7) 


and achieves the boundary conditions (3). The coefficients b k of the expansion (5) are determined by 

the condition that at the data points x = kh the series gives the modified basic data g(kh) or the original 
measurements corrected by a+/3x. Thus: 


n-1 


bi, = 


a=l 


L oot 
g(oh) sin k — 
n 


( 8 ) 


In the absence of noise, the Fourier coefficients bj, b 2 , ..., b m have certain values, but are practically 
zero beyond b m . Then the Fourier synthesis 


g(x) = bj sin - x + b 2 sin — x + ... + b m sin x ( 9 ) 

properly interpolates the function, not only in the data points but at all points of the range. Thus, all 
of the high frequency components of the noise are eliminated. 
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The method of determining the cutoff frequency is different in the present method from that 
proposed by Lanczos. Terms may be added singly by setting k max = -1, -2, -n. In this way, it can 

be determined visually where noise contributions begin to enter. A method for automatic selection of 
k ma x has also been developed in the present computer program, but has been tested for only a few cases. 

In this method the sum 

P = -J=r (b 1 2 + b 2 2 + ... b n 2 ) 1/2 (9) 

\/n-l 


is first formed. The functions 

Gl(m) = (b 1 2 + b 2 2 + ... +b m 2 ) 1/2 (10) 

Vn-1 

are formed for m = 1, 2, ..., n. The value of /I is then subtracted from each of the values of Gl(m), 
resulting in a curve of the type shown by the solid curve in Figure 1. The upper half of this function is 
fitted to a polynomial of degree 3 by the method of least squares and the differences 


Yl(m) = Gl(m) - Y(m) 


(ID 


where Y(m) are the values obtained from the cubic equation, are calculated. The quantities 


A = 


Yl(m) 

Y(m) 


exp. [-bX 2 (m)] 


( 12 ) 


are examined beginning at m = n and proceeding toward smaller m. In equation (12), 

b = 2.303/X max 2 (13) 

Equation (13) was obtained by setting exp(-bX max 2 ) = 0.1. The exponential factor in equation (12) is 

included to damp some possibly large values of A at large m due to small values of Y(m). When A 
exceeds a certain value, in this case set at 0.6, the value of k m for the truncated series is set at 


k max = n - m + 2 (14) 

where n is the total number of coefficients and m is the number of values of A which have been 
examined. The value of 0.6 may require some adjustment as further cases are examined. The original 
function with the noise removed is obtained by reconstituting the data series with the number of terms 
given by k m included. There is some flexibility in choosing the value of A and the extra quantity of 2 

terms in equation (13), since there is a range of terms for which the smoothed function changes very 
little (shown in Figure 2 for a coated 4130 steel sample). 
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EXPERIMENTAL 


Measurements of corrosion potentials (E^orr), corrosion rates (obtained with the polarization 

resistance method) and electrical resistance were made over a period of 45 days for a 4130 steel sample 
coated with Braycote 137 preservative compound (Sample 1) and for a period of 32 days (Sample 2). 
The sample holder employed is shown in Figure 3. The metal specimens were 1.43 cm (9/16 in.) in 
diameter and 0.137 cm thick. The metal samples were smoothed by wet sanding with 220A silicon 
carbide paper and sprayed on one side with a 20 percent by weight solution of Braycote 137 in 1,1,1- 
trichloroethane to a thickness of 0.005 cm (0.002 in.), as measured with a wet thickness gauge. 

Aluminum samples (2219-T87) were prepared by a 15 min immersion in hot alkaline cleaner, 
followed by a 15 min suspension in “Smut-Go” chromate deoxidizer. The samples were then treated 
with Alodine 1200 (conversion coat) for a period of 2 min and sprayed on one side to a measured 
thickness with TT-P-1757 zinc chromate primer. The samples were placed in the sample holder, with 
an area of 1.0 cm^ exposed, and immersed in a test solution consisting of 3.5 percent NaCl buffered at 
pH 5.5 for the entire test period. 

Data for electrical resistance were obtained with the EG&G-PARC Model 356 IR Compensation 
Module in conjunction with the EG&G-PARC Model 3 50 A Corrosion Measurement Console. Data were 
collected daily for the first few days for each sample, after which the frequency of data collection was 
decreased. Values of E^qrr and corrosion current (I^ORR^ were obtained using the polarization 

resistance method where possible, with data being taken on alternate days. The small currents involved 
in the study of coated surfaces disturb the sample surfaces very little, so that repeated measurements 
can be made. 

The EG&G-PARC Model 350A Corrosion Measurement Console was used for collection of 
polarization resistance data at 25°C. Data were collected at 0.5 mV intervals at a scan rate of 0.1 mV/ 
sec. The measurement range for all determinations was -20 to +20 mV with respect to E^qrr, with 

all data being corrected for IR-drop. The data were stored on disk and transferred to a computer for 
calculation of the polarization resistance (R p ), E^qrr, anodic and cathodic Tafel constants and IcORR 
using the program POLCURR [4] . 


RESULTS AND DISCUSSION 


A. Braycote 137 Preservative 

The variation of E^qrr with time for Braycote 137 coated sample 2 is shown in Figure 4. In 
this case, the experimental values of E^qrr, shown at the left side of Figure 4, were obtained through 

least squares refinement of the experimental data by the polarization resistance method using POLCURR. 
The smoothed data are shown at the right. In this case, three terms in the truncated series, as chosen by 
the automatic selection technique, were required for the smoothed data. A peak occurs at about 20 days 
of the exposure period. This is almost exactly the same period for the peak occurrence as obtained by 
Wormwell and Brasher [2] (Fig. 5) for paint-coated steel samples, although their data were the average 
for a large number of samples. The present data are the result of a single measurement for each time 
interval. The variations of corrosion rate and electrical resistance for the same sample are shown in 
Figure 6. A peak in the resistance-time curve occurs at 15 days, while the corrosion rate-time curve 
begins to increase after about 20 days. The E^oRR-time curve thus correlates rather well with the 
curves of Figure 6. 
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The EcQRp-time curves for Braycote 137 sample 1 are shown in Figure 7, with the observed 

values on the left and the smoothed values on the right. In this case, no measurable corrosion current 
was observed until after 27 days, which precluded the use of POLCURR. The measurements were, 
therefore, made manually and resulted in a great deal of data scatter. The smoothed curve at the right 
contained four terms in the truncated series using the automatic selection method. The first peak at 
7 days is not considered significant, and the major peak of the curve occurs at about 31 days. A peak 
occurs at 24 days in the resistance-time curve (Fig. 8), while the corrosion rate becomes measurable at 
about 27 days. The maximum for the smoothed E(-.Qpp-time curve therefore correlates well with the 
curves of Figure 8. 


B. Primer Coated Aluminum 

The EcQpp-time curves for primer coated 2219-T87 aluminum are shown in Figure 9. The 

measurements were carried out over a period of 30 days, during which the curves showed no significant 
variation. This is in contrast to observations for the coated steel samples in the present work, as well as 
the results of Wormwell and Brasher [2]. However, the resistance-time curves and corrosion rate-time 
curves (Fig. 10) for the two samples show a great deal of activity, with a sharp drop in resistance after 
only a few days and a peak in the corrosion rate-time curves at about 16 days. It appears, therefore, 
that E^Qpp measurements cannot be used in the case of aluminum to evaluate corrosion behavior 

unless, possibly, the measurements are extended over a much longer time period. 


CONCLUSIONS 


It is clear from these results that there are important differences in the behavior of coated 
aluminum and steel as far as electrochemical measurements are concerned. The E^Qpp-time curves for 

Braycote 137 coated steel show a maximum after a period of several days, the positions of which corre- 
late well with observations for resistance-time curves and corrosion rate-time curves. Also, since there is 
considerable scatter in the data, a smoothing procedure must be used before proper interpretation of 
the data can be accomplished when measurements are made with a single sample. On the other hand, the 
EcoRR*ti me curves for primer coated aluminum show no significant variations with time over a 30 day 

period, although considerable activity is indicated in the resistance-time and corrosion rate-time curves. 
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Figure 8. Resistance-time and corrosion rate-time curves for Bray cote Sample 1. 
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Figure 9. E^q^r versus time for primer coated 2219-T87 aluminum. 
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APPENDIX 


Listing of the Fortran-77 program used for smoothing the E CO RR-time data. 


PROGRAM main 

ANALYSIS OF FIRST ORDER DIFFERENCE RESULTS 
ORIGINAL DATA MUST BE EQUALLY SPACED. IF THEY ARE NOT, THE 
APPROPRIATE QUESTION IN THE DATA INPUT MUST BE ANSWERED NO. 
THE INCLUDED LAGRANGIAN INTERPOLATING ROUTINE WILL PROVIDE 
EQUAL SPACING WITH THE ADDITION OF ONE DATA POINT. 

KMAX IS DEFINED IN SUBROUTINE LOPAS. THE NUMBER OF 
COEFFICIENTS IN THE TRUNCATED SERIES CAN BE CHOSEN AT WILL 
BY SETTING KMAX EQUAL TO -1,-2, .. . ,-NPTS. AUTOMATIC 
DETERMINATION OF THE NUMBER OF COEFFICIENTS IN THE 
TRUNCATED SERIES IS CARRIED OUT BY A LEAST SQUARES 
FITTING PROCEDURE IF KMAX IS ENTERED AS 1. THIS METHOD 
HAS NOT YET BEEN THOROUGHLY TESTED. 


IMPLICIT REAL*8 ( A-H , 0- Z ) 

CHARACTER* 3 AB, AC 
CHARACTER* 10 IFILEN 
DIMENSION TITLE (18) ,X(100) ,Y(100) 
DIMENSION TX(201) ,TY(201) ,AUX(201) 
COMMON TX.TY.AUX 


SET UP PROBLEM AND INPUT DATA 


15 

5 

25 


30 


75 

125 

20 

200 

310 

320 


215 


220 

330 

C 

C 

C 

C 

C' 

35 


WRITE(6, 15) 

FORMAT(/, IX, ’READ TITLE’,/) 

READ (5, 5) TITLE 
FORMAT (18A4) 

WRITE( 6 , 25) 

FORMAT (/, IX, ’READ LOWER VALUE OF X FOR SMOOTHING’,/) 
READ( 5,20) SL 

WRITE ( 6 , 30 ) „ 

FORMAT (/, IX, ’READ UPPER VALUE OF X FOR SMOOTHING’,/) 
READ ( 5 , 20 ) SU 
WRITE (6 ,75) 

FORMAT (/, IX, ’READ KMAX’,/) 

READ( 5 , 125 ) KMAX 
FORMAT (13) 

FORMAT (F10.0) 

FORMAT ( IX, F8. 4, IX, F9. 5) 

WRITE(6, 310) 

FORMAT (/, IX, ’SAVE SMOOTHED DATA? (YES) (NO) ’ , /) 

FORMATS A3) 

READ ( 5 , 320 ) AB 
IF(AB.EQ. ’YES’ ) THEN 
WRITE(6 , 215) 

FORMAT (/, IX, ’READ FILENAME FOR SMOOTHED DATA’,/) 
READ( 5,220) IFILEN 
END IF 
FORMAT (A10) 

WRI TE (6 330) 

FORMAT (/, IX, ’ARE DATA EQUALLY SPACED? (YES) (NO) ’, /) 
READ( 5 , 320 ) AC 


READ OBSERVED DATA TO SENTINEL. READ ISENT,X,Y:E.G. , 

0,1.0,. 235 CARRIAGE RETURN.' CONTINUE FOR ALL DATA POINTS. 

AFTER ALL DATA HAVE BEEN ENTERED, ENTER 1,0,0 CARRIAGE RETURN. 


WRITE(6 , 35) 

FORMAT (/, IX, ’READ DATA TO SENTINEL’,/) 
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NX=1 

40 READ ( 5 , * ) ISENT,X(NX) , Y(NX) 

IF(ISENT.EQ. 1) GO TO 45 
NX=NX+1 
GO TO 40 
45 NX=NX-1 

NPTS=NX 
C 

C CARRY OUT INTERPOLATION IF DATA POINTS ARE NOT EQUALLY SPACED 

C 

IF(AC.EQ. ’NO’ ) THEN 
DEL=(X(NX)-X(1) ) /FLOAT (NX) 

NPTS=NX+1 
DO 90 1=1 , NPTS 
K=I-1 

TX ( I ) =FLOAT ( K ) *DEL+X ( 1 ) 

90 . TY( I ) =YLAG(TX( I ) , X, Y, 0 , 3 , NX, IEX) 

GO TO 340 
END IF 

DO 345 1=1, NPTS 
TX(I)=X(I) 

345 TY ( I ) =Y ( I ) 

340 OPEN ( 4 , FILE= ’ PRN . LST ’ ) 

WRITE (4, 50) TITLE 
50 FORMAT ( IX , 1 8A4 , / ) 

CALL LOPAS ( SL , SU , NPTS , KMAX ) 

DO 300 K=1 , NPTS 
TY2=TY(K) -AUX(K) 

300 AUX(K) =TY2 

C 

C OUTPUT ORIGINAL AND SMOOTHED DATA 

C 

115 WRITE( 4 , 55 ) 

55 FORMAT(/,7X, ’ORIGINAL DATA’ ,/) 

WRITE( 4 , 60 ) 

60 FORMAT ( 5X, ’X’ , 14X, ’ Y’ , /) 

WRITE (4 ,65)(X(K),Y(K),K=1,NX) 

65 FORMAT (Ell . 5, 3X.E11 . 5) 

WRITE(4 , 70 ) 

70 FORMAT ( / , 7X , ’ SMOOTHED DATA ’ , / ) 

WRITE(4, 60) 

WRITE (4 , 65 ) (TX(K) , AUX(K) , K=1 , NPTS) 

CLOSE (-4) 

IF(AB.EQ. ’YES’ ) THEN 

OPEN ( 2 , FILE=IFILEN , STATUS= ’ NEW’ ) 

WRITE ( 2 , 200 ) (TX( I ) , AUX( I ) , 1=1 , NPTS) 

CLOSE ( 2 , STATUS = ’ KEEP ’ ) 

END IF 
END 
C 

FUNCTION YLAG(XI ,X, Y, IND1 ,N1 , IMAX, IEX) 

IMPLICIT REAL*8(A-H,0-Z) 

c 

C LAGRANGIAN INTERPOLATION 

C 

DIMENSION X( 1 ) , Y( 1 ) 

C 

IND=IND1 

N=N1 

IEX=0 
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IF(N.LE.IMAX) GO TO 10 

N=IMAX 

IEX=N 

10 IF(IND.GT.O) GO TO 40 
DO 20 J=1 , IMAX 
IF(XI-X( J) ) 30,45,20 
20 CONTINUE 

IEX=1 
GO TO 70 
30 IND=J 

40 IF(IND.GT.l) GO TO 50 

IEX=-1 

50 INL=IND- (N+l )/2 

IF(INL.GT.O) GO TO 60 

INL=1 

60 INU=INL+N-1 

IF( INU. LE. IMAX) GO TO 80 

70 ' INL=IMAX-N+1 
INU=IMAX 
80 S=0.0 

P=1.0 

DO 25 J=INL, INU 
P=P*(XI-X(J)) 

D=1 . 0 

DO 15 I=INL, INU 
IF(I.NE.J) GO TO 90 
XD=XI 
GO TO 15 
90 XD=X( J) 

15 D=D*(XD-X( I ) ) 

25 S=S+Y( J)/D 

YLAG=S*P 
35 RETURN 

45 YLAG=Y( J) 

GO TO 35 
END 
C 

SUBROUTINE LOPAS ( SL , SU , NS , KMAX ) 

REMOVES NOISE FROM DATA SET,S,DI BY SMOOTHING IN THE LARGE 
LANCZOS, APPLIED ANALYSIS ( 1956) , P.331-338 

S X-array 

DI Y-array to be smoothed 

NS Number of X,Y points 

BY Noise determined by LOPAS 

COEF Fourier Coefficients 

GX Working array 

KMAX Number of coefficients in truncated series 

>0: KMAX to be found by LOPAS 
<0: KMAX= I ABS ( KMAX ) read in 
SL No smoothing for values X<X(SL) 

SU No smoothing for values X>X(SU) 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION GX( 257 ),V(8),Y(20i),Gl( 100 ) , Y1 ( 201 ) 

DIMENSION S( 201 ) , DI ( 201 ) , BY(201 ) , COEF ( 512) 

COMMON S , DI , BY 

DATA ZRO, TWO, PI/0 . 0D0 , 2 . 0D0 , 3 . 141592654D0/ 

C 
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c 

c 

c 


12 


15 

20 


C 

C 

C 


25 


210 

200 


1 

305 


310 

C 

C 

C 


42 


IF( (SU. LT. SL) .OR. (NS.GT. 257) ) GO TO 100 
EXPAND DATA INTO FOURIER SERIES 
PS=PI/FLOAT(NS) 

SLP= (DI (NS) -DI (1))/(S(NS)-S(1)) 

NN=NS-1 
GX(NS)=ZRO 
GX( 1 ) =ZRO 
COEF( 1 ) =0 
DO 12 N=2 , NN 

GX(N) =DI (N) -DI ( 1 ) -S(N)*SLP 
DO 20 K=1 , NS 
SUM=ZRO 
DO 15 J=2 , NS 

SUM=SUM+GX ( J ) *DSIN ( K* ( J-l ) *PS ) 

COEF (K)= TWO*SUM/FLOAT ( NS ) 

REMOVE HIGH FREQUENCY TERMS 

SUM=ZRO 
DO 25 J=1 , NS 
NJ=NS- J+l 

SUM=SUM+COEF(NJ ) **2 
BETA=DSQRT( SUM/FLOAT ( NS- 1 ) ) 

IF(KMAX. LT.O) THEN 
KMAX=IABS(KMAX) 

T=100 . 0 
GO TO 42 
END IF 

DO 200 1=1, NS 

NI=NS-I+1 

G1 (NI )=0 . 0 

DO 210 J=1 , NI 

TERM=COEF( J)**2 

G1 ( NI ) =G1 ( NI ) +TERM 

G1 ( NI ) =DSQRT ( G1 ( NI ) /FLOAT ( NS-1 ) ) 

G1 ( NI ) =DABS ( G1 ( NI ) -BETA ) 

CALL LSTSQR(NS, S, G1 , T, V) 

DO 305 1=1, NS 
NI=NS-I+1 

yS(*NsI*-S(NI^ j**3 ))+V(3) * ( ^ ( ^ S)- ^ ( ^ I) )**2+V( 4 ) * 

Y1 (NI ) =G1 (NI )-Y(NI ) 

B=2 . 302585/S(NS)**2 
DO 310 1=1, NS 
NI=NS-I+1 

IF(Y(NI) .EQ.0.0) GO TO 310 
A=Y1(NI)/Y(NI)*DEXP(-1.0*B*S(NI)**2) 

IF(A.GT. 0.60) THEN 

KMAX=NI+1 

GO TO 42 

END IF 

CONTINUE 

RECONSTITUTE DATA SERIES 

DO 50 J=1 , NS 
SUM=ZRO 

DO 45 K=1 , KMAX 
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4 5 SUM=SUM+COEF ( K ) *DS I N ( K* ( J - 1 ) *PS ) 

50 BY ( J ) =SUM+DI ( 1 ) +SLP*S ( J ) 

C 

WRITE(4 , 5000 ) BETA, NS , KMAX, SL, SU 
DO 60 K=1 , NS 

IF( (S(K) .LT.SL) .OR. (S(K) .GT.SU ) ) THEN 
BY(K) =0 . 0 
GO TO 60 
END IF 

BY(K)=DI(K)-BY(K) 

60 CONTINUE 

C 

5000 FORMAT(/, IX, ’LOP AS CALLED’ ,/, IX, 'BETA = ’ 

1 , El 1 . 5 , / , IX , ’ NCOF = ’,15, ’ , NTERMS = ’,15,/, IX, ’SL = ’,F5.2, ’ 

1 , SU = ’ , F6 . 2, /) 

C 

WRITE( 4 , 300 ) T 

300 FORMAT(/, IX, ’PERCENT GOODNESS OF FIT = ’,F7.2,/) 

RETURN 

C 

100 WRITE(4, 5500) 

DO 70 K=1 , NS 
70 BY ( K ) = ZRO 

RETURN 

5500 FORMAT ( 1H ,’LOW PASS DATA WINDOW NOT USED’,/) 

END 

C 

SUBROUTINE LSTSQR(NS,S,G1 ,T,V) 

SUBROUTINE FOR FITTING POLYNOMIAL OF DEGREE D 
T=PERCENT GOODNESS OF FIT 

HERE, UPPER HALF OF DATA ARE FITTED BY LEAST SQUARES WITH 
A POLYNOMIAL WITH DEGREE EQUAL TO 3 

IMPLICIT REAL*8(A-H,0-Z) 

INTEGER D,D2 

DIMENSION S(201) ,G1(100) ,A(8,8) ,R(8) ,V(8) ,P(50) ,X(100) ,Y(100) 

C 

NA=NS/2 
NP=NS-NA+1 
D=3 - 

D2=2*D 
N=D+1 

DO 250 1=1, NP 

NI=NS-I+1 

X(I)=S(NS)-S(NI) 

250 Y(I)=G1(NI) 

DO 15 J=2 , D2+1 

P(J)=0.0 

DO 20 K=1 , NP 

20 P(J)=P( J)+X(K)**(J-1) 

15 CONTINUE 

P( 1 ) =NP 
R( 1 ) =0 . 0 
DO 25 J=1 , NP 
25 R( 1 )=R( 1 )+Y( J) 

IF(N.EQ.l) GO TO 30 
DO 35 J=2 , N 
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R( J)=0 . 0 
DO 40 K=1 , NP 

40 R(J)=R(J)+Y(K)*X(K)**(J-1) 

35 CONTINUE 

30 DO 45 J=1 , N 

DO 50 K=1 , N 
50 A(J,K)=P(J +K- 1 ) 

45 CONTINUE 

IF(N.EQ.l) THEN 
V ( 1 )=R( 1 )/A( 1 , 1 ) 

GO TO 110 
END IF 

DO 55 K=1 , N-l 

I=K+1 

L=K 

60 IF ( ( DABS (A(I,K))).GT.( DABS (A(L,K)))) THEN 

L=I 

END IF 

IF(I.LT.N) THEN 
1 = 1+1 
GO TO 60 
END IF 

IF(L.EQ.K) GO TO 210 
DO 65 J=K, N 
Q=A(K, J) 

A(K, J)=A(L, J) 

65 A(L, J)=Q 

Q=R(K) 

R(K)=R(L) 

R(L) =Q 
210 I=K+1 

70 Q=A( I , K)/A(K, K) 

A( I , K)=0 . 0 
DO 75 J=K+1 ,N 

75 A(I,J)=A(I,J)-Q*A(K,J) 

R( I ) =R( I ) -Q*R(K) 

IF(I.LT.N) THEN 
1 = 1+1 
GO TO 70 
END IF 

55 CONTINUE 

V(N)=R(N)/A(N,N) 

DO 80 T=N-1 ,1,-1 
Q=0 . 0 

DO 85 J=I+1,N 
Q=Q+A(I, J)*V(J) 

85 V(I)=(R(I)-Q)/A(I,I) 

80 CONTINUE 

110 Q=0 . 0 

DO 90 J=1 , NP 
90 Q=Q+Y( J) 

FM=Q/FLOAT ( NP ) 

T=0 . 0 
G=0 . 0 

DO 95 J=2 , NP 
Q=0 . 0 

DO 100 K=1,N 

l'OO Q=Q+V(K)*X(J)**(K-1) 

T=T+(Y(J)-Q)**2 
95 G=G+ ( Y ( J ) -FM ) **2 
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IF(G.EQ.O.O) THEN 
T=100 . 0 
GO TO 105 
END IF 

T=100 . 0*DSQRT( 1 . O-T/G) 
10 S RETURN 
END 
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